Data preparation

This document includes code to reproduce all analyses and figures in the paper Nielsen, Dablander, Debnath, Emegor, Ghai, Gwozdz, Hahnel, Hofmann, & Bauer (2024). Perceived plasticity of climate-relevant behaviors and policy support among high- and lower-income individuals.

library(rlang)
library(knitr)
library(Hmisc)
library(dplyr)
library(ggplot2)
library(corrplot)
library(tidyverse)
library(doParallel)
library(BayesFactor)
library(RColorBrewer)
source('helpers.R')

df <- readRDS('Data/BPdata.RDS') %>%
  mutate(
    country = factor(country, c('Denmark', 'United States', 'Nigeria', 'India'))
  )
df$climate_concern <- rowMeans(df[,c('cc_worry', 'cc_importance')], na.rm = TRUE)

actions_order <- c(
  # Curtailment
  'Use less energy for heating and cooling your home',
  'Drive fewer kilometers/miles in your car',
  'Take fewer national and international flights',
  'Eat less red meat',
  'Eat less white meat',
  
  # Investment
  'Improve the insulation of your home',
  'Replace older appliances with newer energy efficient models',
  'Install solar panels at home',
  'Purchase an electric vehicle',
  'Move private investments to climate-friendly financial products'
)

Main Figures

Figure 1: Behavioral plasticity across countries

We look at Behavioral plasticity across countries.

varnames <- c(paste0('curtailment_', seq(5)), paste0('investment_', seq(5)))
map_names <- label(df[varnames])

d <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit()

d <- d %>% 
  group_by(country, name, value) %>% 
  summarize(count = n()) %>% 
  left_join(
    d %>% group_by(country, name) %>% summarize(total_count = n()),
    by = c('country', 'name')
  ) %>% 
  mutate(
    prop = count / total_count
  ) %>% 
  mutate(
    value_name = factor(
      ifelse(value == 1, '1: Very unlikely', ifelse(value == 5, '5: Very likely', value)),
      levels = rev(c('1: Very unlikely', '2', '3', '4', '5: Very likely'))
    ),
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

cols <- c('#ff7f0e', '#A65628', '#029E73', '#56B4E9', '#0072B2')
p_country <- ggplot(d, aes(x = country, y = prop, fill = as.factor(value_name))) +
  geom_bar(stat = 'identity') +
  facet_wrap(~ full_name, ncol = 5) + 
  xlab('') +
  ylab('Proportion') +
  scale_fill_manual(
    values = rev(cols)
    # labels = rev(c('1: Very unlikely', '2', '3', '4', '5: Very likely'))
  ) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    panel.grid = element_blank(),
    legend.title = element_blank(),
    strip.text = element_text(size = 7),
    axis.text.x = element_text(angle = 90, hjust = 1),
    plot.title = element_text(size = 14, hjust = 0.50)
  ) +
  guides(fill = guide_legend(reverse = TRUE))

ggsave('Figures/Figure1.pdf', p_country, width = 9, height = 7)
p_country

Calculate proportions for describing the figure.

library(pander)

d_comp <- d %>% 
  group_by(country, full_name) %>% 
  summarize(
    unlikely = round(sum(prop[value %in% c(1, 2)]), 2),
    likely = round(sum(prop[value %in% c(4, 5)]), 2)
  ) %>% data.frame

pander(d_comp)
country full_name unlikely likely
Denmark Use less energy for heating and cooling your home 0.29 0.45
Denmark Drive fewer kilometers/miles in your car 0.48 0.27
Denmark Take fewer national and international flights 0.43 0.33
Denmark Eat less red meat 0.43 0.38
Denmark Eat less white meat 0.57 0.2
Denmark Improve the insulation of your home 0.42 0.36
Denmark Replace older appliances with newer energy efficient models 0.22 0.53
Denmark Install solar panels at home 0.57 0.25
Denmark Purchase an electric vehicle 0.48 0.31
Denmark Move private investments to climate-friendly financial products 0.46 0.26
United States Use less energy for heating and cooling your home 0.27 0.48
United States Drive fewer kilometers/miles in your car 0.36 0.43
United States Take fewer national and international flights 0.35 0.43
United States Eat less red meat 0.47 0.33
United States Eat less white meat 0.59 0.21
United States Improve the insulation of your home 0.33 0.47
United States Replace older appliances with newer energy efficient models 0.25 0.49
United States Install solar panels at home 0.47 0.37
United States Purchase an electric vehicle 0.5 0.33
United States Move private investments to climate-friendly financial products 0.45 0.27
Nigeria Use less energy for heating and cooling your home 0.22 0.58
Nigeria Drive fewer kilometers/miles in your car 0.35 0.44
Nigeria Take fewer national and international flights 0.3 0.46
Nigeria Eat less red meat 0.28 0.53
Nigeria Eat less white meat 0.44 0.38
Nigeria Improve the insulation of your home 0.14 0.63
Nigeria Replace older appliances with newer energy efficient models 0.09 0.73
Nigeria Install solar panels at home 0.08 0.82
Nigeria Purchase an electric vehicle 0.3 0.52
Nigeria Move private investments to climate-friendly financial products 0.16 0.63
India Use less energy for heating and cooling your home 0.07 0.78
India Drive fewer kilometers/miles in your car 0.14 0.7
India Take fewer national and international flights 0.16 0.64
India Eat less red meat 0.2 0.65
India Eat less white meat 0.21 0.61
India Improve the insulation of your home 0.09 0.74
India Replace older appliances with newer energy efficient models 0.08 0.77
India Install solar panels at home 0.08 0.82
India Purchase an electric vehicle 0.09 0.79
India Move private investments to climate-friendly financial products 0.09 0.73

We use the BayesFactor R package for statistical inference.

Figure 2: Behavioral plasticity and income

We show the relationship between behavioral plasticity and income across items and countries.

df_preds <- do.call('rbind', lapply(varnames, function(varname) get_df_pred_bf(varname, df))) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

country_colors <- c('#DE425B', '#3C3B6E', '#008000', '#FF9933')
p_bp_income <- ggplot(df_preds, aes(y = ypred, x = combined_income, fill = country)) +
  geom_vline(aes(xintercept = 10), linetype = 'dotted', color = 'gray48') +
  geom_point(aes(y = mean_value, x = combined_income, col = country), size = 1, show.legend = FALSE) +
  geom_line(aes(col = country), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40) + 
  xlab('Income group') +
  ylab('Perceived behavioral plasticity') +
  scale_x_continuous(breaks = c(seq(1, 15, 3), 15), limits = c(1, 15)) +
  scale_y_continuous(breaks = seq(1, 5), limits = c(1, 5)) +
  scale_color_manual(values = country_colors) +
  scale_fill_manual(values = country_colors) +
  facet_wrap(~ full_name) +
  theme_minimal() +
  facet_wrap(~ full_name, ncol = 5) +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(size = 8),
    legend.title = element_blank(),
    legend.position = 'top'
  )

ggsave('Figures/Figure2.pdf', p_bp_income, width = 10, height = 6)
p_bp_income

We calculate inclusion Bayes factors in favor of an effect.

df_inclusion_bfs <- get_bf_inclusion_fig2_df(actions_order)
pander(df_inclusion_bfs %>% select(-rscaleFixed, -rscaleCont))
income income_x_country full_name
0.237 0.243 Use less energy for heating and cooling your home
0.396 3.601 Drive fewer kilometers/miles in your car
0.42 5654 Take fewer national and international flights
9.363e+10 0.276 Eat less red meat
2226095 865 Eat less white meat
0.267 4.643 Improve the insulation of your home
7.62 17.21 Replace older appliances with newer energy efficient models
0.154 0.096 Install solar panels at home
3.938e+17 39.3 Purchase an electric vehicle
422.4 0.108 Move private investments to climate-friendly financial products

Figure 3: Behavioral plasticity and domain-matched policy support

We run the relevant three-way interaction models.

df_meat_narm <- df %>%
  filter(!is.na(support_7), !is.na(curtailment_2)) %>% 
  mutate(
    curtailment_2 = curtailment_2 - mean(curtailment_2),
    combined_income = combined_income - mean(combined_income)
  )

df_flights_narm <- df %>%
  filter(!is.na(support_8), !is.na(curtailment_1)) %>% 
  mutate(
    curtailment_1 = curtailment_1 - mean(curtailment_1),
    combined_income = combined_income - mean(combined_income)
  )
  
df_finance_narm <- df %>%
  filter(!is.na(support_6), !is.na(investment_4)) %>% 
  mutate(
    investment_4 = investment_4 - mean(investment_4),
    combined_income = combined_income - mean(combined_income)
  )

rscaleFixed <- sqrt(1/2)
rscaleCont <- sqrt(1/2) / 4

# Lose 422 rows
fit_meat <- generalTestBF(
  support_7 ~ curtailment_2 * combined_income * country,
  data = df_meat_narm, whichModels = 'withmain', 
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Lose 443 rows
fit_flights <- generalTestBF(
  support_8 ~ curtailment_1 * combined_income * country,
  data = df_flights_narm, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Lose 530 rows
fit_finance <- generalTestBF(
  support_6 ~ investment_4 * combined_income * country,
  data = df_finance_narm, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Get the model that predicts the data best
fit_meat_best <- fit_meat[which.max(fit_meat@bayesFactor[, 1])]
fit_flights_best <- fit_flights[which.max(fit_flights@bayesFactor[, 1])]
fit_finance_best <- fit_finance[which.max(fit_finance@bayesFactor[, 1])]

## Used below
countries <- c('Denmark', 'United States', 'Nigeria', 'India')
newdat <- expand.grid(
  country = countries,
  combined_income = seq(15) - mean(seq(15)),
  bp_value = seq(5) - mean(seq(5))
)

map_back <- function(x) {
  as.numeric(factor(x, labels = seq(5)))
}

We look at model predictions using the models with the highest marginal likelihood.

matches <- list(
  c('support_8', 'curtailment_1'),
  c('support_7', 'curtailment_2'),
  c('support_6', 'investment_4')
)

df_preds <- bind_rows(
  get_df_pred_policy(fit_meat_best, 'curtailment_2', 'support_7', df_meat_narm),
  get_df_pred_policy(fit_flights_best, 'curtailment_1', 'support_8', df_flights_narm),
  get_df_pred_policy(fit_finance_best, 'investment_4', 'support_6', df_finance_narm)
) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(map_names)
    ),
    
    country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India'))
  )
matching_colors <- c("#E69F00", "#009E73", "#56B4E9")

p_matching <- ggplot(df_preds, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Perceived behavioral plasticity') +
  ylab('Policy support') +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_y_continuous(breaks = seq(1, 7, 1)) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/Figure3.pdf', p_matching, width = 9, height = 8)
p_matching

library(VGAM)

model_size_prior <- dbetabinom.ab(seq(0, 7), size = 7, shape1 = 1, shape2 = 1)
x <- model_size_prior[1]

# We specify the priors here according to the sequence from the output of the BayesFactor package
# This is if we wanted to use the betabinomial prior, but for the paper we report a uniform model prior
# There's barely any difference
model_prior <- c(
  # Size 0: 1 --> grand mean only
  # Size 1: 3 --> one main effect
  # Size 2: 3 --> two main effects
  # Size 3: 4 --> all main effects or two main effects and one interaction
  # Size 4: 3 --> all main effects and one interactions
  # Size 5: 3 --> all main effects and two two-way interactions
  # Size 6: 1 --> all main effects and all two-way interactions
  # Size 7: 1 --> all main effects and all interactions
  x,
  rep(x / 3, 3),
  rep(x / 3, 3),
  rep(x / 4, 4),
  rep(x / 3, 3),
  rep(x / 3, 3),
  x,
  x
)

set.seed(1)
df_inclusion_bfs3 <- data.frame(rbind(
  get_bf_inclusion_fig3(fit_flights, 'curtailment_1'),
  get_bf_inclusion_fig3(fit_meat, 'curtailment_2'),
  get_bf_inclusion_fig3(fit_finance, 'investment_4')
))

colnames(df_inclusion_bfs3) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
df_inclusion_bfs3$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]

pander(df_inclusion_bfs3)
Table continues below
bp bp_country bp_income bp_income_country
2.011e+88 651.9 10.94 0.003961
9.906e+149 123.9 4.204 0.01475
1.36e+101 6800 6124 0.004789
full_name
Take fewer national and international flights
Eat less red meat
Move private investments to climate-friendly financial products

Not so different with uniform prior, which we report in the main paper.

set.seed(1)
df_inclusion_bfs3 <- data.frame(rbind(
  get_bf_inclusion_fig3(fit_flights, 'curtailment_1', prior = 'uniform'),
  get_bf_inclusion_fig3(fit_meat, 'curtailment_2', prior = 'uniform'),
  get_bf_inclusion_fig3(fit_finance, 'investment_4', prior = 'uniform')
))
colnames(df_inclusion_bfs3) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
df_inclusion_bfs3$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]

pander(df_inclusion_bfs3)
Table continues below
bp bp_country bp_income bp_income_country
1.661e+88 681.3 11.43 0.003961
8.183e+149 129.5 4.394 0.01475
1.124e+101 7107 6400 0.004789
full_name
Take fewer national and international flights
Eat less red meat
Move private investments to climate-friendly financial products

Extended Data Figures

Figure ED1: Already did the action or cannot do it

We visualize the proportion of people who already did a particular action or cannot do it.

varnames_peb <- paste0(varnames, '_peb')
varnames_cannot <- paste0(varnames, '_cannot')
map_names_cannot <- map_names
names(map_names_cannot) <- varnames_cannot

df_cannot <- df %>%
  select(country, all_of(varnames_cannot)) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    prop = mean(value == 1), # Only 0 / 1 exists
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_cannot[name]),
      levels = wrap_label_vectorized(actions_order)
    ),
    prop = ifelse(grepl('curtailment', name), NA, prop),
    prop_se = ifelse(grepl('curtailment', name), NA, prop_se)
  )

varnames_peb <- paste0(varnames, '_peb')
map_names_peb <- map_names
names(map_names_peb) <- varnames_peb

df_peb <- df %>%
  select(country, all_of(varnames_peb)) %>%
  pivot_longer(-country) %>% 
  filter(str_starts(name, 'curtailment')) %>% 
  mutate(value = as.numeric(as.character(value))) %>% 
  group_by(country, name) %>% 
  summarize(
    prop = mean(value == -1), 
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  bind_rows(
    df %>%
      select(country, all_of(varnames_peb)) %>%
      pivot_longer(-country) %>% 
      mutate(value = as.numeric(as.character(value))) %>% 
      filter(str_starts(name, 'investment')) %>% 
      group_by(country, name) %>% 
      summarize(
        prop = mean(value == 1), 
        prop_se = sqrt(prop * (1 - prop) / n())
      )
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_peb[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

df_both <- bind_rows(
  df_cannot %>% mutate(type = 'Cannot do'),
  df_peb %>% mutate(type = 'Already did / Never do it')
) %>% 
  mutate(
    
    # We change the labels
    full_name = case_when(
      full_name == 'Use less energy for\nheating and cooling your\nhome' ~ 'Do not heat or cool their home',
      full_name == 'Drive fewer kilometers/miles in\nyour car' ~ 'Do not own a car',
      full_name == 'Eat less red meat' ~ 'Do not eat red meat',
      full_name == 'Eat less white meat' ~ 'Do not eat white meat',
      full_name == 'Take fewer national and\ninternational flights' ~ 'Never fly',
      TRUE ~ full_name
    ),
    
    full_name = factor(
      full_name,
      levels = c(
        rev(unname(wrap_label_vectorized(actions_order)))[seq(5)],
        'Do not eat white meat',
        'Do not eat red meat',
        'Never fly',
        'Do not own a car',
        'Do not heat or cool their home'
      )
    )
  )

pwidth <- 0.96
p_both <- ggplot(df_both, aes(x = prop, y = full_name, fill = country, col = country)) +
  geom_bar(stat = 'identity', position = position_dodge(width = pwidth)) +
    geom_point(
    aes(x = prop, y = full_name), position = position_dodge(width = pwidth),
    size = 2, show.legend = FALSE, color = 'black'
  ) +
  geom_errorbar(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.40, linewidth = 1, show.legend = FALSE, color = 'black'
  ) +
  geom_point(
    aes(x = prop, y = full_name), position = position_dodge(width = pwidth),
    size = 1, show.legend = FALSE
  ) +
  geom_errorbar(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.30, linewidth = 0.30,
    show.legend = FALSE
  ) +
  facet_wrap(~ type) +
  ylab('') +
  xlab('Proportion') +
  scale_fill_manual(values = country_colors) +
  scale_color_manual(values = country_colors) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 10),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureED1.pdf', p_both, width = 9, height = 8)
p_both

Figure ED2: Already did and behavioral plasticity

Visualize relationship between exclusion rate and mean plasticity rating for behaviors across countries.

df_mean <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    se_bp = sd_bp / sqrt(n())
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

df_mean_exclusion <- df_mean %>% 
  left_join(
    df_both %>%
      filter(type == 'Already did / Never do it') %>% 
      mutate(
        name = str_extract(name, "^[^_]+_[^_]+")
      ) %>% 
      select(country, name, prop, prop_se)
    , by = c('country', 'name')
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = rev(wrap_label_vectorized(actions_order))
    )
  )
  
cols <- brewer.pal(10, 'Set3')
p_mean_exclusion <- ggplot(df_mean_exclusion, aes(x = prop, y = mean_bp, col = full_name)) +
  geom_point(size = 2) +
  geom_errorbar(
    aes(ymin = mean_bp - 1.96 * se_bp, ymax = mean_bp + 1.96 * se_bp),
    width = 0.01
  ) +
  # Horizontal error bars
  geom_errorbarh(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    height = 0.1
  ) +
  facet_wrap(~ country) +
  labs(
    x = 'Proportion: Already do / never do it',
    y = 'Mean perceived behavioral plasticity'
  ) +
  scale_color_manual(name = '', values = cols) +
  scale_x_continuous(breaks = seq(0, 0.3, 0.05)) +
  scale_y_continuous(breaks = seq(1, 5, 1), limits = c(1, 5)) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    strip.text = element_text(size = 14),
    plot.title = element_text(hjust = 0.50, size = 14),
    axis.text.y = element_text(size = 10),
    legend.key.spacing = unit(0.10, 'cm')
  ) +
  guides(color = guide_legend(nrow = 5, reverse = TRUE))

ggsave('Figures/FigureED2.pdf', p_mean_exclusion, width = 7, height = 7)
p_mean_exclusion

Figure ED3: Already did the action across income groups

df_prep <- df %>%
  select(country, combined_income, varnames_peb) %>%
  pivot_longer(cols = -c(country, combined_income)) %>% 
  mutate(value = as.numeric(as.character(value))) %>% 
  mutate(
    income_bracket = ifelse(
      combined_income <= 5, 'low',
      ifelse(combined_income < 10, 'medium', 'high')
    )
  )

df_selection <- df_prep %>% 
  filter(str_starts(name, 'curtailment')) %>% 
  group_by(country, income_bracket, name) %>% 
  summarize(
    prop = mean(value == -1), 
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  bind_rows(
    df_prep %>%
      filter(str_starts(name, 'investment')) %>% 
      group_by(country, income_bracket, name) %>% 
      summarize(
        prop = mean(value == 1), 
        prop_se = sqrt(prop * (1 - prop) / n())
      )
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_peb[name]),
      levels = wrap_label_vectorized(actions_order)
    ),
    
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    ),
    
    # We change the labels
    full_name = case_when(
      full_name == 'Use less energy for\nheating and cooling your\nhome' ~ 'Do not heat or cool their home',
      full_name == 'Drive fewer kilometers/miles in\nyour car' ~ 'Do not own a car',
      full_name == 'Eat less red meat' ~ 'Do not eat red meat',
      full_name == 'Eat less white meat' ~ 'Do not eat white meat',
      full_name == 'Take fewer national and\ninternational flights' ~ 'Never fly',
      TRUE ~ full_name
    ),
    
    full_name = factor(
      full_name,
      levels = c(
        'Do not heat or cool their home',
        'Do not own a car',
        'Never fly',
        'Do not eat red meat',
        'Do not eat white meat',
        rev(unname(wrap_label_vectorized(actions_order)))[seq(5)]
      )
    )
  )

p_selection <- ggplot(
    df_selection,
    aes(y = prop, x = country, fill = income_bracket)
  ) +
  geom_bar(stat = 'identity', position = position_dodge(width = pwidth)) +
  geom_point(
    aes(x = country, y = prop), position = position_dodge(width = pwidth),
    size = 2, show.legend = FALSE, color = 'black'
  ) +
  geom_errorbar(
    aes(ymin = prop - 1.96 * prop_se, ymax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.40, linewidth = 1, show.legend = FALSE, color = 'black'
  ) +
  geom_point(
    aes(x = country, y = prop, col = income_bracket), position = position_dodge(width = pwidth),
    size = 1, show.legend = FALSE
  ) +
  geom_errorbar(
    aes(ymin = prop - 1.96 * prop_se, ymax = prop + 1.96 * prop_se, col = income_bracket),
    position = position_dodge(width = pwidth),
    width = 0.30, linewidth = 0.30,
    show.legend = FALSE
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('Proportion') +
  xlab('') +
  scale_fill_manual(values = matching_colors) +
  scale_color_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 7),
    axis.text.x = element_text(angle = 90),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureEx3.pdf', p_selection, width = 9, height = 7)
p_selection

Let’s run a proportion test to assess the evidence for differences.

df_bfs <- c()
for (country in unique(df_prep$country)) {
  for (name in unique(df_prep$name)) {
    bf <- test_prop(df, country, name)
    df_bfs <- rbind(df_bfs, cbind(bf, country, name))
  }
}

df_bfs <- data.frame(df_bfs) %>% 
  mutate(
    name = gsub('_peb$', '', name),
    log_bf = as.numeric(bf),
    full_name = map_names[name]
  )

pander(df_bfs %>% select(country, full_name, log_bf))
country full_name log_bf
Denmark Take fewer national and international flights 45.66
Denmark Eat less red meat -5.937
Denmark Eat less white meat -6.544
Denmark Drive fewer kilometers/miles in your car 24.79
Denmark Use less energy for heating and cooling your home -5.283
Denmark Purchase an electric vehicle 25.19
Denmark Replace older appliances with newer energy efficient models 3.275
Denmark Improve the insulation of your home 22.6
Denmark Move private investments to climate-friendly financial products -1.516
Denmark Install solar panels at home -0.9868
United States Take fewer national and international flights 52.32
United States Eat less red meat -5.661
United States Eat less white meat -6.602
United States Drive fewer kilometers/miles in your car 17.95
United States Use less energy for heating and cooling your home -6.944
United States Purchase an electric vehicle -0.01785
United States Replace older appliances with newer energy efficient models 0.4505
United States Improve the insulation of your home 6.468
United States Move private investments to climate-friendly financial products -5.966
United States Install solar panels at home 0.8944
India Take fewer national and international flights -1.688
India Eat less red meat -4.436
India Eat less white meat -4.312
India Drive fewer kilometers/miles in your car -1.241
India Use less energy for heating and cooling your home -5.221
India Purchase an electric vehicle 1.757
India Replace older appliances with newer energy efficient models -4.226
India Improve the insulation of your home -3.277
India Move private investments to climate-friendly financial products -4.846
India Install solar panels at home -4.477
Nigeria Take fewer national and international flights 7.376
Nigeria Eat less red meat -5.844
Nigeria Eat less white meat 3.438
Nigeria Drive fewer kilometers/miles in your car 8.633
Nigeria Use less energy for heating and cooling your home -1.016
Nigeria Purchase an electric vehicle -4.526
Nigeria Replace older appliances with newer energy efficient models 25.45
Nigeria Improve the insulation of your home -4.727
Nigeria Move private investments to climate-friendly financial products -5.887
Nigeria Install solar panels at home 30.04

Supplementary Figures

Figure S1: Mean behavioral plasticity across countries

Here we visualize mean behavioral plasticity across all countries.

df_mean_sd_beh <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  mutate(curtailment = as.numeric(grepl('curtailment', name))) %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    n = n(),
    se_bp = sd_bp / sqrt(n)
  ) %>% 
  # Compute the mean of the mean /sd behavioral plasticity of behaviors
  group_by(country) %>% 
  summarize(
    mean_bp = mean(mean_bp),
    sd_bp = mean(sd_bp)
  )

pander(df_mean_sd_beh)
country mean_bp sd_bp
Denmark 2.777 1.365
United States 2.902 1.461
Nigeria 3.563 1.292
India 4.024 1.149
df_mean <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    n = n(),
    se_bp = sd_bp / sqrt(n),
    se_sd_bp = sd_bp / sqrt(2 * n)
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

pwidth <- 0.91
p_mean <- ggplot(df_mean, aes(x = mean_bp, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = mean_bp - 1.96 * se_bp, xmax = mean_bp + 1.96 * se_bp),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Mean perceived behavioral plasticity') +
  scale_x_continuous(breaks = seq(1, 5, 1), limits = c(1, 5)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 7),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS1.pdf', p_mean, width = 9, height = 7)
p_mean

Figure S1: Standard deviation of behavioral plasticity across countries

pwidth <- 0.91
p_sd <- ggplot(df_mean, aes(x = sd_bp, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = sd_bp - 1.96 * se_sd_bp, xmax = sd_bp + 1.96 * se_sd_bp),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Standard deviation of perceived behavioral plasticity') +
  scale_x_continuous(breaks = seq(0, 2, 0.50), limits = c(0, 2)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 7),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS2.pdf', p_sd, width = 9, height = 7)
p_sd

Figure S3: Mean policy support across countries

Here we visualize mean policy support across all countries.

support_vars <- colnames(df %>% select(starts_with('support_')))
support_labels <- label(df %>% select(starts_with('support_')))

df_mean_ps <- df %>%
  select(country, starts_with('support_')) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_ps = mean(value),
    sd_ps = sd(value),
    n = n(),
    se_ps = sd_ps / sqrt(n)
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(support_labels[name]),
      levels = wrap_label_vectorized(support_labels[name])
    )
  )

pwidth <- 0.91
country_colors <- c('#DE425B', '#3C3B6E', '#008000', '#FF9933')
p_mean_ps <- ggplot(df_mean_ps, aes(x = mean_ps, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = mean_ps - 1.96 * se_ps, xmax = mean_ps + 1.96 * se_ps),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Mean policy support') +
  scale_x_continuous(breaks = seq(1, 7, 1), limits = c(1, 7)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 6),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS3.pdf', p_mean_ps, width = 9, height = 8)
p_mean_ps

Figure S4: Differences in high and low income groups

Here we visualize the differences in mean behavioral plasticity between the top 10% income earners and the rest.

set.seed(1)
df_diffs <- do.call(
  'rbind', lapply(varnames, function(varname) get_top10_difference(varname, df))
) %>% 
  mutate(
    country = factor(country, levels = rev(countries)),
    full_name = factor(
      map_names[name],
      levels = rev(actions_order)
    )
  )
cols <- brewer.pal(10, 'Set3')
p_estimates <- ggplot(df_diffs, aes(x = mean_diff, y = country, color = full_name)) +
  geom_vline(aes(xintercept = 0), linetype = 'dotted', color = 'gray48') +
  geom_errorbar(
    aes(xmin = mean_diff_lo, xmax = mean_diff_hi),
    position = position_dodge(width = 0.70), width = 0.80, size = 1.5
  ) +
  geom_point(position = position_dodge(width = 0.70), size = 3) + 
  scale_x_continuous(breaks = round(seq(-0.80, 1.2, 0.2), 3), limits = c(-0.90, 1.3)) +
  theme_minimal() +
  scale_color_manual(name = '', values = cols) +
  labs(
    x = 'Mean difference in perceived behavioral plasticity',
    y = ''
  ) +
  theme(
    legend.position = 'top',
    strip.text = element_text(size = 14),
    plot.title = element_text(hjust = 0.50, size = 14),
    axis.text.y = element_text(size = 10),
    legend.key.spacing = unit(0.06, 'cm')
  ) +
  guides(color = guide_legend(nrow = 5, reverse = TRUE))

ggsave('Figures/FigureS4.pdf', p_estimates, width = 8, height = 10)
p_estimates

Figures S5 & S6: Sensitivity analysis I

We assess how sensitive the results are to changes in the prior width of the income effect and the interaction of income and country.

scale_lo <- 0.05
scale_hi <- 0.30
scales_cont <- seq(scale_lo, scale_hi, 0.0125)

df_inclusion_bfs_sens <- do.call('rbind', lapply(scales_cont, function(scale) {
  get_bf_inclusion_fig2_df(actions_order, rscaleCont = scale)
}))
end <- Sys.time()

df_inclusion_bfs_sens <- df_inclusion_bfs_sens %>% 
    mutate(
    full_name = factor(
      wrap_label_vectorized(full_name),
      levels = wrap_label_vectorized(actions_order)
    )
  )
library(scales)

# Custom function to apply scientific notation for numbers >= 10
custom_label <- function(x) {
  ifelse(x > 100, scientific_format()(x), format(x, scientific = FALSE))
}

df_inclusion_bfs_sens2 <- df_inclusion_bfs_sens %>% 
  mutate(income = ifelse(income <= 1, 1 / income, income))

p_sens_fig2_income <- ggplot(df_inclusion_bfs_sens, aes(x = rscaleCont, y = (income))) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = custom_label) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  facet_wrap(~ full_name, ncol = 5, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS5.pdf',
  p_sens_fig2_income,
  width = 10, height = 5
)
p_sens_fig2_income

p_sens_fig2_income_country <- ggplot(
  df_inclusion_bfs_sens, aes(x = rscaleCont, y = (income_x_country))
) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = custom_label) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  facet_wrap(~ full_name, ncol = 5, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS6.pdf',
  p_sens_fig2_income_country,
  width = 10, height = 5
)
p_sens_fig2_income_country

Figure S7: Sensitivity analysis II

We run sensitivity analyses for the results reported in Figure 3.

if (!file.exists('Data/df_sensitivity_matched.csv')) {
  
  rscaleFixed <- sqrt(1/2)
  scales_cont <- seq(0.05, 0.30, 0.0125)
  
  df_inclusion_bfs <- do.call('rbind', lapply(scales_cont, function(scale) {
    
    set.seed(1)
    fit_meat <- generalTestBF(
      support_7 ~ curtailment_2 * combined_income * country,
      data = df_meat_narm, whichModels = 'withmain', 
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    set.seed(1)
    fit_flights <- generalTestBF(
      support_8 ~ curtailment_1 * combined_income * country,
      data = df_flights_narm, whichModels = 'withmain',
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    set.seed(1)
    fit_finance <- generalTestBF(
      support_6 ~ investment_4 * combined_income * country,
      data = df_finance_narm, whichModels = 'withmain',
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    df_inclusion_bfs <- data.frame(rbind(
      get_bf_inclusion_fig3(fit_flights, 'curtailment_1', prior = 'uniform'),
      get_bf_inclusion_fig3(fit_meat, 'curtailment_2', prior = 'uniform'),
      get_bf_inclusion_fig3(fit_finance, 'investment_4', prior = 'uniform')
    ))
    
    colnames(df_inclusion_bfs) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
    df_inclusion_bfs$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]
    df_inclusion_bfs$rscaleCont <- scale
    df_inclusion_bfs
    
  }))
  
  df_inclusion_bfs_sens2 <- df_inclusion_bfs %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(full_name),
        levels = wrap_label_vectorized(actions_order)
      )
    )
  
  df_inclusion_bfs_sens2_long <- df_inclusion_bfs_sens2 %>%
    pivot_longer(
      cols = c(bp, bp_country, bp_income, bp_income_country),
      names_to = 'type',
      values_to = 'value'
    ) %>% 
    mutate(
      type = case_when(
        type == 'bp' ~ 'BP',
        type == 'bp_country' ~ 'BP x Country',
        type == 'bp_income' ~ 'BP x Income',
        type == 'bp_income_country' ~ 'BP x Income x Country'
      )
    )
} else {
  df_inclusion_bfs_sens2_long <- read.csv('Data/df_sensitivity_matched.csv') %>% 
    mutate(
      full_name = factor(
        full_name, levels = c(
          'Take fewer national and\ninternational flights', 'Eat less red meat',
          'Move private investments to\nclimate-friendly financial products'
        )
      )
    )
}
p_sens_fig3 <- ggplot(
  df_inclusion_bfs_sens2_long, aes(x = rscaleCont, y = value)
) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = label_scientific()) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  # facet_nested(~ type + full_name, scales = 'free_y') +
  facet_wrap(~ type + full_name, ncol = 3, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS7.pdf',
  p_sens_fig3,
  width = 8, height = 8
)
p_sens_fig3

Figures S8 & S9: More loosely domain-matched analyses

We estimate the relationship between more loosely domain-matched policies and behavioral plasticity items.

get_best <- function(fit_all) fit_all[which.max(fit_all@bayesFactor[, 1])]

matches_loose <- list(
  
  # Increasing energy efficiency requirements for buildings & insulating your home
  c('support_5', 'investment_3'),
  
  # Increase price of electricity during peak hours & replace older with more energy efficient appliances
  c('support_3', 'investment_2'),
  
  # increasing subsidies for renewable energy projects & install solar panels at home
  c('support_4', 'investment_5'),
  
  # Ban sale of diesel & purchasing an EV
  c('support_10', 'investment_1'),
  
  # Requires banks to disclose & move investments
  c('support_6', 'investment_4'),
  
  # Increase subsidies for food products with low GHGs & eating less red meat
  c('support_12', 'curtailment_2'),
  
  # Increase subsidies for food products with low GHGs & eating less red meat
  c('support_12', 'curtailment_3'),
  
  # Expand public transport & drive fewer kilometres with the car
  c('support_2', 'curtailment_4'),
  
  # Strengthen requirements for energy efficiency & use less energy for heating and cooling
  c('support_5', 'curtailment_5')
)

#' Prepare data.frame for loosely domain-matched policy analysis
prepare_df <- function(df, bp_name, policy_name) {
  bp_name <- ensym(bp_name)
  policy_name <- ensym(policy_name)
  
  df %>% 
    filter(!is.na(!!bp_name), !is.na(!!policy_name)) %>% 
    mutate(
      !!bp_name := !!bp_name - mean(!!bp_name, na.rm = TRUE),
      combined_income = combined_income - mean(combined_income, na.rm = TRUE)
    )
}

if (!file.exists('Data/df_preds_loose.csv')) {
  
  fit_cars <- get_bf_fit_loose(df, 'support_10', 'investment_1')
  fit_electricity <- get_bf_fit_loose(df, 'support_3', 'investment_2')
  fit_banks <- get_bf_fit_loose(df, 'support_6', 'investment_4')
  fit_replace <- get_bf_fit_loose(df, 'support_5', 'investment_3')
  fit_renewables <- get_bf_fit_loose(df, 'support_4', 'investment_5')
  
  fit_food_red <- get_bf_fit_loose(df, 'support_12', 'curtailment_2')
  fit_food_white <- get_bf_fit_loose(df, 'support_12', 'curtailment_3')
  fit_public <- get_bf_fit_loose(df, 'support_2', 'curtailment_4')
  fit_efficiency <- get_bf_fit_loose(df, 'support_5', 'curtailment_5')
  
  df_preds_loose <- bind_rows(
    get_df_pred_policy(
      get_best(fit_cars), 'investment_1', 'support_10',
      prepare_df(df, 'investment_1', 'support_10')
    ) %>% mutate(type = 'cars'),
    
    get_df_pred_policy(
      get_best(fit_electricity), 'investment_2', 'support_3',
      prepare_df(df, 'investment_2', 'support_3')
    ) %>% mutate(type = 'electricity'),
    
    get_df_pred_policy(
      get_best(fit_banks), 'investment_4', 'support_6',
      prepare_df(df, 'investment_4', 'support_6')
    ) %>% mutate(type = 'banks'),
    
    get_df_pred_policy(
      get_best(fit_replace), 'investment_3', 'support_5',
      prepare_df(df, 'investment_3', 'support_5')
    ) %>% mutate(type = 'replace'),
    
    get_df_pred_policy(
      get_best(fit_renewables), 'investment_5', 'support_4',
      prepare_df(df, 'investment_5', 'support_4')
    ) %>% mutate(type = 'renewables'),
    
    get_df_pred_policy(
      get_best(fit_food_red), 'curtailment_2', 'support_12',
      prepare_df(df, 'curtailment_2', 'support_12')
    ) %>% mutate(type = 'red_meat'),
    
    get_df_pred_policy(
      get_best(fit_food_white), 'curtailment_3', 'support_12',
      prepare_df(df, 'curtailment_3', 'support_12')
    ) %>% mutate(type = 'white_meat'),
    
    get_df_pred_policy(
      get_best(fit_public), 'curtailment_4', 'support_2',
      prepare_df(df, 'curtailment_4', 'support_2')
    ) %>% mutate(type = 'public'),
    
    get_df_pred_policy(
      get_best(fit_efficiency), 'curtailment_5', 'support_5',
      prepare_df(df, 'curtailment_5', 'support_5')
    ) %>% mutate(type = 'public')
  ) %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(map_names[name]),
        levels = wrap_label_vectorized(actions_order)
      ),
      country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India'))
    )
} else {
  df_preds_loose <- read.csv('Data/df_preds_loose.csv') %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(map_names[name]),
        levels = wrap_label_vectorized(actions_order)
      ),
      country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India')),
      income_bracket = factor(income_bracket, levels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)'))
    )
}

df_preds_curtailment_loose <- df_preds_loose %>% 
  filter(str_starts(name, 'curtailment'))

df_preds_investment_loose <- df_preds_loose %>% 
  filter(str_starts(name, 'investment'))
matching_colors <- c("#E69F00", "#009E73", "#56B4E9")

p_investment_loose <- ggplot(
  df_preds_investment_loose, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Perceived behavioral plasticity') +
  ylab('Policy support') +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  scale_y_continuous(breaks = seq(1, 7), limits = c(1, 7)) +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

p_curtailment_loose <- ggplot(
  df_preds_curtailment_loose, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Behavioral plasticity') +
  ylab('Policy support') +
  scale_y_continuous(breaks = seq(1, 7), limits = c(1, 7)) +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/FigureS8.pdf', p_curtailment_loose, width = 9, height = 10)
p_curtailment_loose

ggsave('Figures/FigureS9.pdf', p_investment_loose, width = 9, height = 10)
p_investment_loose

Figure S10: Mean behavioral plasticy and mean policy support

We run the model for mean support across all policies and mean of all behavioral plasticity items.

df_avg <- df %>% 
  rowwise() %>% 
  mutate(
    support_avg = mean(c_across(starts_with('support')), na.rm = TRUE),
    bp_value = mean(c_across(c(paste0('curtailment_', seq(5)), paste0('investment_', seq(5)))), na.rm = TRUE)
  ) %>% data.frame

fit_all <- generalTestBF(
  support_avg ~ bp_value * combined_income * country,
  data = df_avg, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

fit_all_best <- fit_all[which.max(fit_all@bayesFactor[, 1])]

df_preds_avg <- get_df_pred_policy(
  fit_all_best, 'bp_value', 'support_avg', df_avg, type = 'averaged',
  join = FALSE, use_empirical = FALSE
) %>% 
  mutate(
    country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India')),
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    )
  )

df_emp_avg <- df_avg %>% 
  mutate(
    income_bracket = ifelse(
      combined_income <= 5, 'low',
      ifelse(combined_income < 10, 'medium', 'high')
    ),
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    )
  )

df_both <- df_emp_avg %>% 
  left_join(df_preds_avg, by = c('country', 'bp_value', 'income_bracket'))

p_mean_matched <- ggplot(df_both, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = support_avg, col = income_bracket),
    size = 1, alpha = 0.30, width = 0.10, show.legend = FALSE
  ) +
  geom_line(data = df_preds_avg, aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(data = df_preds_avg, aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  scale_y_continuous(breaks = seq(1, 7, 1), limits = c(0.9, 7.2)) +
  xlab('Mean perceived behavioral plasticity') +
  ylab('Mean policy support') +
  facet_wrap(~ country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/FigureS10.pdf', p_mean_matched, width = 9, height = 4)
p_mean_matched

Supplementary Figures 11 - 14: Correlation plots

We calculate Kendall’s correlation coefficients.

support_vars <- colnames(df %>% select(starts_with('support_')))
support_labels <- label(df %>% select(starts_with('support_')))

df_cors <- df %>% select(all_of(varnames), starts_with('support_'), country)
labels <- c(wrap_label_vectorized(map_names), wrap_label_vectorized(support_labels))

res_india <- get_cors_all(df_cors, labels, 'India')
res_denmark <- get_cors_all(df_cors, labels, 'Denmark')
res_nigeria <- get_cors_all(df_cors, labels, 'Nigeria')
res_usa <- get_cors_all(df_cors, labels, 'United States')

Denmark

pdf('Figures/FigureS11.pdf', width = 10, height = 10)
get_corrplot_all(res_denmark)
dev.off()
## quartz_off_screen 
##                 2
get_corrplot_all(res_denmark)

United States

pdf('Figures/FigureS12.pdf', width = 10, height = 10)
get_corrplot_all(res_usa)
dev.off()
## quartz_off_screen 
##                 2
get_corrplot_all(res_usa)

Nigeria

pdf('Figures/FigureS13.pdf', width = 10, height = 10)
get_corrplot_all(res_nigeria)
dev.off()
## quartz_off_screen 
##                 2
get_corrplot_all(res_nigeria)

India

pdf('Figures/FigureS14.pdf', width = 10, height = 10)
get_corrplot_all(res_india)
dev.off()
## quartz_off_screen 
##                 2
get_corrplot_all(res_india)